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Introduction 


Numerical methods for viscous transonic aerodynamics have made enormous 
progress in recent years, but further improvement is still needed, especially with 
regard to the accuracy of the prediction. This improvement will mainly result from 
validation of codes based on data from appropriate transonic experiments where 
e.g. boundary conditions to be used in the calculation have to be measured and 
where also shear stress data or velocity profiles are available. Because there are 

"°< e "° U0 . h SUCh ex P erimen,s H is also necessary to compare numerical results 
with data from experiments which don't meet the requirements of code validation. 
Even In these cases a lot of useful information can be gained, not onlv for thp 
theoretician but dlso for the experimentalist. 

In the following a comparison is presented between the computed results of the 
flow about the CAST10 airfoil and the pressure distributions and force coefficents 
from experiments in the adaptive TCT [1], 

Description of the Method 

For the discretization of the complete Navier-Stokes equations a finite volume 
Runge-Kutta time-stepping scheme based on a cell centered formulation is used 
[2]. Since the fluxes across the cell faces are averaged from neighbouring cells 
which is equivalent to central differencing artificial diffusion terms are needed. 
These damping terms are the usual blend of second and fourth order differences 
a wei 9 hting function is employed. This weighting function restricts the 
artificial diffusion in the viscous near-wall regions or In the wake, where the 
unweighted formulation would otherwise result in a predominance of non-physical 
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The set of ordinary differential equations resulting from the finite volume discreti- 
zation is integrated in time using a linearized four-stage Runge Kutta scheme. 
Local time stepping is employed to accelerate the convergence to steady-state 
solutions and in order to save further computation time the artificial and physical 
diffusion terms are updated only once per time step reducing the execution time 
by more than 50 per cent. 


To have a well-posed problem a set of appropriate boundary conditions is needed. 
At the airfoil surface the no-slip and the adiabatic wall condition is used. The 
pressure i$_derived from the assumption of zero pressure gradient normal to the 
wall which is justified for the very small step sizes normal to the wall used in 
Navier-Stokes calculations for turbulent flows. At the far field boundary we use 

one-dimensional Riemann invariants normal to the boundary in order to obtain 
boundary conditions. 


Since the two-dimensional flow solver used here is derived from the three-dimen- 
sional one described in [2], it allows for a block-structured approach, i.e. the 
computational domain can be divided in a number of subdomains called blocks. 
The advantage of this approach is that it is very flexible regarding the handlinq 
of complex geometries (e.g. with multiply connected domains). The flexibility is 
partly due to the segmentation of the block faces, which allows for the use of dif- 
ferent types of boundary conditions and different neighbouring blocks at each 



block face. This feature was used in the present calculations to model the finite 
thickness of the CAST10 trailing edge by using a C-type mesh past the airfoil with 
an extra grid in the gap behind the trailing edge. For more Information on the 
block structure refer to [2]. 

To simulate the turbulence the well-known algebraic model of Baldwin and Lomax 
[3] is used. 



Typical grid In the vicinity of the airfoil 

For the flow about the airfoil the mesh consists of 260x80 cells in chordwise and 
wall-normal direction, respectively and the mesh extends about 10 root chords 
away from the airfoil. There are 200 cells on the wing surface and the first step size 
normal to the wing surface is chosen such that it is equivalent to z*x 2 in the first 
cell of the wall, in order to resolve the laminar sublayer of the turbulent flow. 

The C-type grid is generated algebraically using a code of Sobieczky [4] resulting 
in a mesh with a gap of the thickness of the trailing edge. This gap is closed by a 
suitable interpolation which yields a smooth distribution of the stepsize in the 
direction normal to the wake. 

The block structure in the present calculations includes three blocks; the first and 
second corresponding to the C-lype mesh in the vicinity of the airfoil and to the 
interpolated mesh in the gap, respectively. In these two blocks the complete 
Navier-Stokes equations are solved, whereas only the Euler equations are solved 
in the third block for the outer part of the C-mesh. The block boundary of the later 

block is about 25 per cent of chord away from the airfoil and the wake respec- 
lively. ’ 


Details of the mesh near the leading and trailing edges 

The figures show the strong clustering of the grid lines near the wall and give an 
idea of the grid in the gap behind the trailing edge. In the present calculation the 
mesh within the gap behind the trailing edge is fairly small, i.e. there are only ten 
cells over the height of the trailing edge. 
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Comparison of experimental and computed lift and drag 

M = .73 Re = 10 million 

The table shows a comparison of the lift and drag coefficients from experiments in 
he TCT LU with free transition (the transition takes place somewhere near the 
leading edge) and calculations where the transition point Is prescribed. The tran- 
sition location is given as upper/lower chord position. 

As can be seen from the table the upper surface location of transition (.07 to .2 
chord) has almost no influence on the coefficients, but the lift increases when the 
transition location is shifted downstream on the lower surface ( to .15 and .2 chord). 

The calculation for a = 3.0 is slightly unsteady, i.e. the residuals stay on a certain 
level in the separation zone behind the shock. 

For details see the following pages. 


a 

exp. /cal. 

transition 

lift coeff. 

drag coeff. 

1.0 

exp. 

free 

.616 

.010 


cal. 

.07/. 07 

.677 

.0152 


cal. 

.10/10 

.677 

.0152 


cal. 

.15/. 15 

.685 

.0149 


cal. 

.07/20 

.695 

.0150 


cal. 

.20/07 

.677 

.0150 


cal. 

.10/.10 

(modif.turb.) 

.637 

.0159 

.9 

cal. 

.10/.10 

.677 

.0150 

3.0 

exp. 

free 

.895 

.0450 

— 

cal. 

.10/. 10 

.82.90 

.033-.039 




Iso-Mach contours for M= .73 Re= 10 million « = 1.00 

Calculation: transition at .101.10 chord 

This figure shows contour lines of the Mach number in the vicinity of the airfoil 
(A= .05) with the supersonic regime set off in gray. Also shown are the block 
boundaries of the computation. Since some isolines cross the boundary between 
the viscous inner and the inviscid outer block without any disturbance it is obvi- 
ous that the block concept has no influence on the solution quality. 
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Cp-distribution at M= .73 Re= 10 million * = 1.00 

1. Calculation: transition at .101.10 (or .071.07) chord 

The symbols show the experimental data from TCT [1](point 203) and the solid line 
shows the computational result. The horizontal line with the crosses at its ends 
indicates the critical Cp value. Fixing the transition at 10 or 7 per cent chord in the 
calculation gives the same pressure distribution. 

The results compare quite well for the major part of the surface except near the 

trailing edge, but there is a discrepancy in the shock location.of about 10 percent 
chord. 

To see the influence of the chosen transition location on the results the followina 
variations are made. y 
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Cp-distrlbutlon M= .73 Re= 10 million « = 1.00 

2. Calculation: transition at .071.20 chord 

The pressure distribution shows almost no variation compared to previous calcu- 
lation. 

3. Calculation: transition at .201.07 chord 

The pressure at the beginning of the pressure plateau is now even closer to the 
experimental one, but the shock has moved a little bit further downstream. The 
solution behind the shock and at the lower surface are the same as before. 
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Cp-distribution M— .73 Re= 10 million a = 1.00 

4. Calculation: transition at .101.10 chord and a = .9 


A decrease In the angle of attack results In less agreement regarding the plateau 
UDstream nf P . res * Ure a * the lower surface. The shock has now moved 

hlahPr M«rh b V f3r ^ ! ma " eXlent To ma,ch ,he pla,eau Pressure again a 
higher Mach number would be necessary, but this further variation was not tried 
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Cp-distributlon M= .73 Re= 10 million * - 1.00 

5. Calculation: transition at .101.10 chord and variation of turbulence model near 
shocks 

A modification in the Baldwin/Lomax turbulence model resulting in a local 
increase of the eddy viscosity is made. The shock is moved upstream almost into 
the right position but the plateau pressure and the pressure at the lower surface 
show now larger discrepancies similar to the result with the smaller angle of 
attack It is obvious from the movement of the shock that there is a strong influence 
of the modelling on the turbulent shock boundary layer interaction. 
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Iso-Mach contours for M= .73 Re= 10 million a = 1.00 

Calculation: transition at .101.10 chord 
Modified turbulence model 

This figure shows the Mach number distribution with the modified turbulence 
model. In comparison to the previous results one realizes the upstream movement 
of the shock due to this modification. Away from the supersonic region the two 

solutions look very similar as one could expect already from the pressure distrib- 
utions. 


Cp-distributlon M= .73 Re= 10 million a = 3.00 

Calculation: transition at .10/. 10 chord 

For this higher angle of attack case ( point 207 of TCT data[1]) the calculation 
indicates a slightly unsteady solution in the separated region on the upper surface 
near the trailing edge. Again the results compare quite well for the major part of 
the surface, but we find again the discrepancy in the shock location. In this case 
with separated flow the aforementioned modification of the turbulence model 
shows almost no influence on the solution, maybe because the modification is only 
local at the shock and does not extend over the whole separated region. 
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Iso-Mach contours for M= .73 Re= 10 million « = 3.00 

Calculation: transition at .10/. 10 chord 

This figure shows the Mach number distribution for the higher angle of attack case, 
where the shock has moved downstream. Due to the separation behind the shock 
the boundary layer has thickened considerably as can be seen in comparison to 
the Mach number distribution for a = 1 . 
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Comparison of experimental and computed lift and drag at M .765 

The table shows a comparison of the force coefficients from experiments with fixed 
and free (for the higher Re numbers) transition to those from calculations where 
the transition location is always at .07/.07. Results are presented for different 
angles of attack and different Re numbers. 

For more details regarding the calculated results see the following pages. 



Experiment 

Calculation 

a 

lift / drag 

lift / drag 

Re= 4 million 

0 

.378 / .012 

.53 / .0169 

.5 

- 

.604 / .0216 

1. 

.57 / .020 

.53- .65 / .019-.028 

Re = 10 million 

0 

.45 / .011 

.53 / ,0155 

.5 

.588 / .015 

.575 / .0185 

1. 

.623 / .024 

,58-.63 / .019-0.25 

Re = 40 million 

0 

.538 / .012 

.41 -.43 / . 010- 011 



Cp-distribution M= .765 Re= 4 million <x = 0.00 

Transition at .07/ .07 chord (calculation) and .06/. 06 chord (experiment) 

For this higher Mach number case (point 39 of TCT data [1]) the numerical result 
compares not so well with the experimental data, not even at the lower surface. 
The computation exhibits a pressure plateau with an expansion peak In front of the 
shock, whereas the experiment shows a double structure of weaker shocks. The 
calculated pressure distribution results in a higher lift (and drag) and In a higher 
trailing edge pressure. 

According to previous experience the difference in the transition location between 
calculation and experiment is estimated to have practically no influence. 
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Cp-dlstribution M= .765 Re= 10 million * = 0.00 

Transition at .071.07 chord (calculation) and .061.06 chord (experiment) 

Increasing the Reynolds number results in almost no change In the computed 
pressure distribution. The experimental pressure (point 79 in TCT data [1]) is now 
slightly higher on the lower surface and slightly lower on the upper surface thus 
producing a higher lift due to the reduced decambering by the thinner boundary 
layers. Now the pressure distributions compare better except for the region of the 
shocks where we find again the double shock in the experimental data. 

It is far from clear why the calculation at the lower Reynolds number doesn't show 
the decambering effect found in the experiment. 
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Cp-distribution M= .765 Re= 40 million * = 0.00 

Transition at .071.07 chord (calculation) and .06/. 06 chord (experiment) 

A further increase in Reynolds number changes the situation in the experiment 
(point 284), l.e. the double shocks merge and form a stronger single shock down- 
stream. Again the pressure is increased at the lower and decreased at the upper 
surface. Although the qualitative result of the computation compares now better to 
the experimental pressure distribution, the quantitative result is much worse. This 
is due to the poor resolution of the very thin boundary layer. With a better resol- 
ution,. however, the computed shock position is again found downstream of the 
experimental one whereas the plateau pressure and the pressure at the lower 
surface are recovered. As for the lower Mach number cases we assume that this 
effect is at least partly due to the turbulent shock boundary layer interaction which 
is not correctly modelled by the Baldwin-Lomax turbulence model. 



78 



Cp-dlstribution M= .765 Re= 4 million a = 1.00 

Transition at .071.07 chord (calculation) and .061.06 chord (experiment) 

For a higher angle of attack and the low Reynolds number (point 40 of TCT data 
[1]) we find only a single shock in the experimental data and, as is seen in most 
of the other cases, the numerical result compares quite well except at the shock. 
But this depends on the picked iteration cycle where the results are plotted, as 
will be discussed on the next page. 
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Cp-distribution M= .765 Re= 10 million a = 1.00 

Transition at .071.07 chord (calculation) and .061.06 chord (experiment) 

Increasing the Reynolds number again (point 339 of TCT data [1]) has little influ- 
ence in the experimental data; only the shock is shifted downstream a little bit. 
One would expect the same for the computed results keeping in mind the results 
for zero incidence. But as is seen from the variation of lift and drag in the pre- 
ceding table there is an unsteadiness in the numerical results at this angle of 
attack, i.e. the solutions do not converge to a steady state. The result shown here 
was obviously taken at a moment where the lift in the calculation was low, whereas 
the result on the previous page corresponds to a situation where the lift was high. 

Since the numerical method uses local time stepping as an acceleration technique 
the unstead iness cannot be interpreted in a physically meaningful way, although 
It indicates that a time accurate calculation at this angle of attack would yield an 
unsteady flow behavior, too. 
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Concluding Remarks 

Results of the simulation of the viscous flow past the CAST10 airfoil have been 
shown for different flow conditions. Since the experiments provide only surface 
pressures and force coefficents the comparison to the numerical results relies on 
these. 

Good agreement of the results is found for the lower Mach number cases except 
for the shock position. As numerical experiments indicate, this seems to be due to 
the turbulent shock boundary layer interaction which is not correctly modelled by 
the algebraic turbulence model employed. 

For the lower Mach number case the influence of the transition location has been 
investigated, too. Changing the transition location at the lower surface has much 
more influence on the pressure distribution than changing it on the upper side. 

For the higher Mach number case the double shock structure found in the exper- 
iment for the lower Reynolds numbers was not reproduced by the numerical sol- 
utions. The reason for this is unknown though it may be due to the turbulence 
modelling. For the higher Reynolds number a better resolution of the boundary 
layer is needed in the computation in order to recover the experimental pressure 
plateau; but then the shock position is still found downstream of the experimental 
one. 
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